# First replication file for Weidmann & Salehyan, Violence and Ethnic Segregation: A Computational Model Applied to Baghdad, ISQ (2013) 57, 52-64
# This file replicates the results presented in Figures 3 and 4
# It requires the fast_sweep_* files to be presented in the working directory
# Java code for the model is available upon request

# Nils B. Weidmann (nils.weidmann@uni-konstanz.de)

setwd("/Users/nilsw/OldProjects/Baghdad/R")

# read data files
seed1 <- read.table("fast_sweep_1.txt", header=T)
seed2 <- read.table("fast_sweep_2.txt", header=T)
seed3 <- read.table("fast_sweep_3.txt", header=T)
res.all <- rbind(seed1, seed2, seed3)
res.all <- subset(res.all, propInsurgents==0.02 & classThreshold==0.7) # use results for proportion of insurgents = 0.02 and a classification threshold of 0.7 for "homogenous" neighborhoods

res.0 <- subset(res.all, bestfitAt>5) # drop runs that never really improved above the starting conditions (best fit later than step 5)
res.best.0 <- subset(res.0, bestfitviolence<=quantile(res.0$bestfitviolence, 0.5) & bestfitethnic<=quantile(res.0$bestfitethnic, 0.5)) # good runs are those where both violence and ethnic distribution can be approximated better than average. This set is \Theta' in the paper.

# Figure 3: proportion of coethnics
pdf("propCoethnics.pdf", width=7, height=7)
plot(density(res.all$am_propCoethnics), lwd=3, col="grey", ylim=c(0, 0.2), xlab=expression(paste(alpha[1], " (proportion co-ethnics)")), main="", cex.lab=1.5)
lines(density(res.best.0$am_propCoethnics), lty=1, lwd=3)
dev.off()

# Figure 4a: experienced attacks
pdf("expAttacks.pdf", width=7, height=7)
plot(density(res.all$mm_expAttacks), lwd=3, col="grey", ylim=c(0, 0.15), xlab=expression(paste(beta[1], " (experienced violence)")), cex.lab=1.5, main="")
lines(density(res.best.0$mm_expAttacks), lty=1, lwd=3)
dev.off()

# Figure 4b: spatially lagged attacks
pdf("expAttacksNb.pdf", width=7, height=7)
plot(density(res.all$mm_expAttacksNb), lwd=3, col="grey", ylim=c(0, 0.07), xlab=expression(paste(beta[2], " (violence, spatial lag)")), cex.lab=1.5, main="")
lines(density(res.best.0$mm_expAttacksNb), lty=1, lwd=3)
dev.off()




